home *** CD-ROM | disk | FTP | other *** search
/ Mac OS 9 Serial Number Archive / SN Archive 2023.11.04.toast / BSNG / SDK / BSNG SDK 2.5 / Libraries / UltraLibU / UltraLibU.c < prev    next >
Encoding:
C/C++ Source or Header  |  1998-02-05  |  15.5 KB  |  553 lines  |  [TEXT/CWIE]
  1. /*
  2.  * $Workfile:: UltraLibU.c                                                    $
  3.  * $Revision:: 1                                                              $
  4.  *
  5.  * $Author:: Buck Rogers                                                      $
  6.  * $Modtime:: 30.09.1997 17:50 Uhr                                            $
  7.  *
  8.  * $History:: UltraLibU.c                                                     $
  9.  * 
  10.  * *****************  Version 1  *****************
  11.  * User: Buck Rogers  Date: 30.09.1997   Time: 18:31 Uhr
  12.  * Created in $/BSNG/Plugins/BSNG SDK/Libraries/UltraLibU
  13.  * Adding subproject 'BSNG' to '$/'
  14.  *
  15.  * $NoKeywords::                                                              $
  16.  */
  17.  
  18.  
  19. /**************************   ULTRALibU(niversal)   **********************
  20.  
  21. This is a Macintosh/Motorola implementation of the Ultra pseudo-random
  22. number generator --  the state-of-the-art as of 11/92.
  23.  
  24. References --
  25.  
  26.     Creator:    Prof. George Marsaglia and Co-workers
  27.                 Dept. of Statistics
  28.                 Florida State University
  29.                 
  30.     Code by:    Dr. Michael P. McLaughlin
  31.                 MITRE Corporation
  32.                 MS W337
  33.                 1820 Dolley Madison Blvd.
  34.                 McLean, VA  22102
  35.                 
  36.                 E-mail: mpmcl@mitre.org
  37.                 
  38.     Please acknowledge the work of others by leaving this header intact.
  39.     
  40. **************************************************************************
  41.     
  42.             ---------------  General Remarks  ---------------
  43.             
  44. The code was written in a Metrowerks™ environment and is compatible with
  45.     either 68K (68020 and higher) or PowerPC environments.
  46. Assembly language was used for one critical low-level function.
  47. Things to watch out for include the following:
  48.  
  49.     sizeof(double) is compiler-option dependent and is here unspecified.
  50.     Longs and floats are 4 bytes.
  51.     Shorts are 2 bytes.
  52.     
  53.     Functions returning bytes or bits are pre-cast to short.
  54.     
  55.     Registers not saved:
  56.     
  57.         68K:    D0-D2, A0-A1
  58.  
  59.         PPC:    r3-r10
  60.     
  61.     WARNING!  Be sure to compile and run UltraUTest in your desired environment.
  62.     Ascertain that your output is identical to UltraLibUTest.out.org, supplied
  63.     with this package.
  64.      
  65.                 ---------------  Timings  ---------------
  66.  
  67. The execution times for the functions supplied are somewhat variable and
  68. are increased by factors such as cache overflow and program segmentation.
  69. The timings below are average values for a standard Macintosh Quadra 800™
  70. and a PowerMac 9500/120™.  The compiler was Metrowerks CW8™ with full
  71. optimization for speed and both the calling module and UltraLibU in the same
  72. segment.  68040 "Native Floating Point" was used in the 68K case.  To
  73. facilitate comparisons, the calling overhead (to a void function with no
  74. parameters) has been subtracted out.  Your timings may differ.
  75.  
  76.               Function                                Microseconds per Call
  77.                                                               68K            PPC   
  78.  
  79.             Ultra_long32                                1.5            0.27
  80.             Ultra_long31                                 1.5            0.27
  81.         
  82.             Ultra_short16                                  1.0            0.24
  83.             Ultra_short15                                  1.0            0.23
  84.  
  85.             Ultra_short8                                 0.77            0.21
  86.             Ultra_short8u                                 0.72            0.21
  87.             Ultra_short7                                  0.75            0.21            
  88.  
  89.             Ultra_short1                                  0.65            0.23
  90.  
  91.             Ultra_uni                                     4.7            0.57
  92.             Ultra_vni                                     4.8            0.61
  93.  
  94.             Ultra_duni                                     6.4            0.91
  95.             Ultra_dvni                                     6.4            0.90
  96.  
  97.             Ultra_norm                                  24.0            2.74
  98.             Ultra_expo                                  28.1            2.53
  99.  
  100. *************************************************************************/
  101.  
  102. #include <math.h>
  103.  
  104. struct {                    /* to restart Ultra from a known beginning */
  105.         float                Ugauss;
  106.         unsigned long    UFSR[37],USWB[37],Ubrw,Useed1,Useed2;
  107.         long                Ubits;
  108.         short                Ubyt,Ubit;
  109.         char                *Uptr;
  110.     } Ultra_Remember;    /* 324 bytes */
  111.  
  112. double            Ultra_2n63,Ultra_2n31,Ultra_2n7;
  113. float                Ultra_gauss;        /* remaining Ultra_norm variate */
  114. unsigned long    Ultra_FSR[37],        /* final random numbers */
  115.                     Ultra_SWB[37],        /* subtract-with-borrow output */
  116.                     Ultra_brw,            /* either borrow(68K) or ~borrow(PPC) */
  117.                     Ultra_seed1,        /* seeds MUST be initialized with */
  118.                     Ultra_seed2;        /*    values > 0 */
  119. long                Ultra_bits;            /* bits for Ultra_short1 */
  120. short                Ultra_byt,            /* # bytes left in Ultra_FSR[37] */
  121.                     Ultra_bit;            /* # bits left in Ultra_bits */
  122. char                *Ultra_ptr;            /* running pointer to Ultra_FSR[] */
  123.  
  124. /*    Ultra_Refill() is the core of Ultra (see Marsaglia and Zaman, 1991,
  125.     "A New Class of Random Number Generators", Annals of Applied Probability,
  126.     vol. 1(3), 426-480). It refills Ultra_SWB[37] via a subtract-with-borrow
  127.     PRNG then superimposes a multiplicative congruential PRNG to produce
  128.     Ultra_FSR[37] which supplies all of the Ultra-random bytes.
  129.     
  130.     The contents of Ultra_FSR[] are EXTREMELY random, even at the bit level,
  131.     and ALL bits appear to be random, not just the most significant bits.
  132.     The period of this compound generator exceeds 10^356. */
  133.  
  134. /** Public Prototypes ***************************************************/
  135.  
  136. long Ultra_long32();
  137. long Ultra_long31();
  138.  
  139. short Ultra_short16();
  140. short Ultra_short15();
  141.  
  142. short Ultra_short8();
  143. short Ultra_short8u();
  144. short Ultra_short7();
  145.  
  146. short Ultra_short1();
  147.  
  148. float Ultra_uni();
  149. float Ultra_vni();
  150.  
  151. double Ultra_duni();
  152. double Ultra_dvni();
  153.  
  154. float Ultra_norm(float mu,float sigma);
  155. float Ultra_expo(float mu);
  156.  
  157. Boolean Ultra_Init();
  158.  
  159. void Ultra_SaveStart();
  160. void Ultra_RecallStart();
  161.  
  162. /** Private Prototype ***************************************************/
  163.  
  164. void Ultra_Refill();
  165.  
  166. /************************************************************************/
  167.  
  168. #ifdef __POWERPC__
  169. #define    ULTRABRW    0xFFFFFFFF
  170. asm void Ultra_Refill()
  171. {
  172.         lwz        r3,Ultra_brw        /* fetch global addresses from TOC */
  173.         lwz        r6,Ultra_SWB
  174.         lwz        r4,0(r3)                /* ~borrow */
  175.         la            r7,48(r6)            /* &Ultra_SWB[12] */
  176.         
  177.         sub        r5,r5,r5                /* clear entire word */
  178.         mr            r8,r5                    /* counter */
  179.         li            r5,1
  180.         sraw        r4,r4,r5                /* restore XER|CA */
  181.         
  182.         li            r8,24
  183.         mtctr        r8
  184.         la            r4,-4(r6)
  185. UR1:    lwzu        r9,4(r7)
  186.         lwz        r10,4(r4)
  187.         subfe        r9,r10,r9            /* r9 -= r10 */
  188.         stwu        r9,4(r4)
  189.         bdnz+        UR1
  190.         mr            r7,r6                    /* &Ultra_SWB */
  191.         li            r8,13
  192.         mtctr        r8
  193.         la            r7,-4(r6)
  194. UR2:    lwzu        r9,4(r7)
  195.         lwz        r10,4(r4)
  196.         subfe        r9,r10,r9            /* r9 -= r10 */
  197.         stwu        r9,4(r4)
  198.         bdnz+        UR2
  199.  
  200.         lwz        r4,0(r3)                /* ~borrow again */
  201.         addme        r4,r5                    /* r5 = 1 */
  202.         neg        r4,r4
  203.         stw        r4,0(r3)                /* new ~borrow */
  204.  
  205.         la            r6,-4(r6)            /* &SWB[-1] */
  206.         lwz        r7,Ultra_FSR
  207.         lwz        r5,Ultra_ptr
  208.         lwz        r4,Ultra_seed1
  209.         stw        r7,0(r5)                /* reset running pointer to FSR */
  210.         la            r7,-4(r7)            /* overlay congruential PRNG */
  211.         lis        r10,1                    /* r10 = 69069 */
  212.         addi        r10,r10,3533
  213.         lwz        r5,0(r4)                /* Ultra_seed1 */
  214.         li            r8,37
  215.         mtctr        r8
  216. UR3:    lwzu        r9,4(r6)                /* SWB */
  217.         mullw        r5,r5,r10            /* Ultra_seed1 *= 69069 */
  218.         xor        r9,r9,r5
  219.         stwu        r9,4(r7)                
  220.         bdnz+        UR3
  221.         stw        r5,0(r4)                /* save Ultra_seed1 for next time */
  222.         lwz        r7,Ultra_byt
  223.         li            r5,148                /* 4*37 bytes */
  224.         sth        r5,0(r7)                /* reinitialize */
  225.         blr
  226. }
  227. #else
  228. #define    ULTRABRW    0x00000000
  229. asm void Ultra_Refill()
  230. {
  231.         machine    68020
  232.  
  233.         MOVE.L    A2,-(SP)                /* not scratch */
  234.         LEA        Ultra_SWB,A2        /* Ultra_SWB[0] */
  235.         LEA        52(A2),A1            /* Ultra_SWB[13] */
  236.         MOVEQ        #0,D0                    /* restore extend bit */
  237.         SUB.L        Ultra_brw,D0
  238.         MOVEQ        #23,D2                /* 24 of these */
  239. UR1:    MOVE.L    (A1)+,D0
  240.         MOVE.L    (A2),D1
  241.         SUBX.L    D1,D0
  242.         MOVE.L    D0,(A2)+
  243.         DBRA        D2,UR1
  244.         LEA        Ultra_SWB,A1
  245.         MOVEQ        #12,D2                /* 13 of these */
  246. UR2:    MOVE.L    (A1)+,D0
  247.         MOVE.L    (A2),D1
  248.         SUBX.L    D1,D0                    /* subtract-with-borrow */
  249.         MOVE.L    D0,(A2)+
  250.         DBRA        D2,UR2
  251.         MOVEQ        #0,D0
  252.         MOVE.L    D0,D1
  253.         ADDX        D1,D0                    /* get borrow bit */
  254.         MOVE.L    D0,Ultra_brw        /*    and save it */
  255.         LEA        Ultra_SWB,A1
  256.         LEA        Ultra_FSR,A2
  257.         MOVE.L    A2,Ultra_ptr        /* reinitialize running pointer */
  258.         MOVE.L    Ultra_seed1,D0
  259.         MOVE.L    #69069,D1            /* overlay congruential PRNG */
  260.         MOVEQ        #36,D2                /* 37 of these */
  261. UR3:    MOVE.L    (A1)+,(A2)
  262.         MULU.L    D1,D0
  263.         EOR.L        D0,(A2)+
  264.         DBRA        D2,UR3
  265.         MOVE.L    D0,Ultra_seed1        /* save global for next time */
  266.         MOVE        #148,Ultra_byt        /* 4*37 bytes left */
  267.         MOVE.L    (SP)+,A2                /* restore */
  268.         RTS
  269. }
  270. #endif
  271.  
  272. /*    Ultra_long32() returns a four-byte integer, U[-2147483648,2147483647]. It may,
  273.     of course, be cast to unsigned long. */
  274.  
  275. long Ultra_long32()
  276. {
  277.     register long    result;
  278.     
  279.     if (Ultra_byt < 4) Ultra_Refill();
  280.     result = *((long *) Ultra_ptr);
  281.     Ultra_ptr += 4; Ultra_byt -= 4;
  282.     return(result);
  283. }
  284.  
  285. /*    Ultra_long31() returns a four-byte integer, U[0,2147483647]. */
  286.  
  287. long Ultra_long31()
  288. {
  289.     register long    result;
  290.     
  291.     if (Ultra_byt < 4) Ultra_Refill();
  292.     result = *((long *) Ultra_ptr);
  293.     Ultra_ptr += 4; Ultra_byt -= 4;
  294.     return(result & 0x7FFFFFFF);
  295. }
  296.  
  297. /*    Ultra_short16() returns a two-byte integer, U[-32768,32767]. */
  298.  
  299. short Ultra_short16()
  300. {
  301.     register short    result;
  302.     
  303.     if (Ultra_byt < 2) Ultra_Refill();
  304.     result = *((short *) Ultra_ptr);
  305.     Ultra_ptr += 2; Ultra_byt -= 2;
  306.     return(result);
  307. }
  308.  
  309. /*    Ultra_short15() returns a two-byte integer, U[0,32767]. */
  310.  
  311. short Ultra_short15()
  312. {
  313.     register short    result;
  314.     
  315.     if (Ultra_byt < 2) Ultra_Refill();
  316.     result = *((short *) Ultra_ptr);
  317.     Ultra_ptr += 2; Ultra_byt -= 2;
  318.     return(result & 0x7FFF);
  319. }
  320.  
  321. /*    Ultra_short8() returns a two-byte integer, U[-128,127].  It gets a random
  322.     byte and casts it to short.  This operation extends the sign bit.
  323.     Consequently, you may NOT cast this function to unsigned short/long (see 
  324.     Ultra_short8u() below). */
  325.  
  326. short Ultra_short8()
  327. {
  328.     register short    result;
  329.     
  330.     if (Ultra_byt < 1) Ultra_Refill();
  331.     result = (short) *Ultra_ptr;
  332.     Ultra_ptr += 1; Ultra_byt -= 1;
  333.     return(result);
  334. }
  335.  
  336. /*    Ultra_short8u() returns a two-byte integer, U[0,255].  It proceeds as in
  337.     Ultra_short8() but clears the high byte instead of extending the sign bit. */
  338.  
  339. short Ultra_short8u()
  340. {
  341.     register short    result;
  342.     
  343.     if (Ultra_byt < 1) Ultra_Refill();
  344.     result = (short) *Ultra_ptr;
  345.     Ultra_ptr += 1; Ultra_byt -= 1;
  346.     return(result & 0xFF);
  347. }
  348.  
  349. /*    Ultra_short7() returns a two-byte integer, U[0,127]. */
  350.  
  351. short Ultra_short7()
  352. {
  353.     register short    result;
  354.     
  355.     if (Ultra_byt < 1) Ultra_Refill();
  356.     result = (short) (*Ultra_ptr & 0x7F);
  357.     Ultra_ptr += 1; Ultra_byt -= 1;
  358.     return(result);
  359. }
  360.  
  361. /*    Ultra_short1() returns a two-byte integer, U[0,1], i.e., a Boolean variate.
  362.     It calls Ultra_long32() and returns the bits one at a time. */
  363.  
  364. short Ultra_short1()
  365. {
  366.     register short    result;
  367.     
  368.     if (Ultra_bit <= 0) {
  369.         Ultra_bits = Ultra_long32();
  370.         Ultra_bit = 32;
  371.     }
  372.     result = (short)(Ultra_bits < 0);            /* return sign bit */
  373.     Ultra_bits += Ultra_bits;                        /* shift left by one */
  374.     Ultra_bit--;
  375.     return(result);
  376. }
  377.  
  378. /*    Ultra_uni() returns a four-byte float, U(0,1), with >= 25 bits of precision.  This
  379.     precision is achieved by continually testing the leading byte, b, of the mantissa.
  380.     If b == 0, it is replaced with a new random byte and the decimal point readjusted.
  381.     This procedure simultaneously ensures that Ultra_uni() never returns zero. */
  382.  
  383. float Ultra_uni()
  384. {
  385.     register double    fac = Ultra_2n31;
  386.     register long        along;
  387.     register short        extra;
  388.     
  389.     along = Ultra_long31();
  390.     if (along >= 0x01000000) return((float)(fac*along));
  391.     for (extra=0;!extra;) {                        /* will not be an infinite loop */
  392.          extra = Ultra_short7();
  393.          fac *= Ultra_2n7;
  394.     }
  395.     along |= (((long)extra) << 24);
  396.     return((float)(fac*along));
  397. }
  398.  
  399. /*    Ultra_vni() returns a four-byte float, U(-1,1), with the same features as
  400.     described above for Ultra_uni(). */
  401.     
  402. float Ultra_vni()
  403. {
  404.     register double    fac = Ultra_2n31;
  405.     register long        along,limit = 0x01000000;
  406.     register short        extra;
  407.     
  408.     if ((along = Ultra_long32()) >= limit)
  409.         return((float)(fac*along));
  410.     else if (-along >= limit)
  411.         return((float)(fac*along));
  412.     for (extra=0;!extra;) {
  413.          extra = Ultra_short7();
  414.          fac *= Ultra_2n7;
  415.     }
  416.     if (along >= 0) {
  417.         along |= (((long)extra) << 24);
  418.         return((float)(fac*along));
  419.     }
  420.     along = -along;
  421.     along |= (((long)extra) << 24);
  422.     return((float)(-fac*along));
  423. }
  424.  
  425. /*    Ultra_duni() and Ultra_dvni() return double-precision U[0,1) and U(-1,1).  In
  426.     both cases, zero IS a remote possibility.  These functions are intended for
  427.     those occasions when seven significant figures are not enough.  If you need
  428.     TYPE double, but not double PRECISION, then it is much faster to use
  429.     Ultra_uni() or Ultra_vni() and cast -- implicitly or explicitly. */
  430.  
  431. double Ultra_duni()
  432. {
  433.     return(Ultra_long31()*Ultra_2n31 + ((unsigned long) Ultra_long32())*Ultra_2n63);
  434. }
  435.  
  436. double Ultra_dvni()
  437. {
  438.     return(Ultra_long32()*Ultra_2n31 + ((unsigned long) Ultra_long32())*Ultra_2n63);
  439. }
  440.  
  441. /*    Ultra_norm() returns a four-byte float, N(mu,sigma), where mu and sigma are the
  442.     mean and standard deviation, resp., of the parent population.  The normal variates
  443.     returned are exact, not approximate.  Ultra_norm() uses Ultra_vni() so there is no
  444.     possibility of a result exactly equal to mu.  Note that mu and sigma must also be
  445.     floats, not doubles. */
  446.  
  447. float Ultra_norm(float mu, float sigma)
  448. {
  449.     register float    fac,rsq,v1,v2;
  450.  
  451.     if ((v1 = Ultra_gauss) != 0.0) {        /* Is there one left? */
  452.         Ultra_gauss = 0.0;
  453.         return(sigma*v1 + mu);
  454.     }
  455.     do {
  456.         v1 = Ultra_vni();
  457.         v2 = Ultra_vni();
  458.         rsq = v1*v1 + v2*v2;
  459.     } while (rsq >= 1.0);
  460.     fac = sqrt(-2.0*log(rsq)/rsq);
  461.     Ultra_gauss = fac*v2;                    /* Save the first N(0,1) */
  462.     return(sigma*fac*v1 + mu);                /* and return the second. */
  463. }
  464.  
  465. /*    Ultra_expo() returns a four-byte float, Exponential(mu).  The parameter, mu, is
  466.     both the mean and standard deviation of the parent population.  It must be a
  467.     float greater than zero. */
  468.  
  469. float Ultra_expo(float mu)
  470. {
  471.     return(-mu*log(Ultra_uni()));
  472. }
  473.  
  474. /*    Ultra_SaveStart() and Ultra_RecallStart() save and restore, respectively, the
  475.     global variables of UltraLibU, from Ultra_gauss through Ultra_ptr.  If you
  476.     think it may be necessary to recall a sequence of random numbers exactly, FIRST
  477.     call Ultra_SaveStart().  To recover the sequence later, call Ultra_RecallStart().
  478.     If you wish to terminate the program and still recover the same random numbers,
  479.     you must save the array Ultra_Remember[] to a file and read it back upon restart.
  480.     For convenience, Ultra_Remember[] is declared in UltraU.h as an array of unsigned
  481.     longs even though this is not really the case. */
  482.  
  483. void Ultra_SaveStart()
  484. {
  485.     register    i;
  486.  
  487.     Ultra_Remember.Ugauss = Ultra_gauss;
  488.     Ultra_Remember.Ubits = Ultra_bits;
  489.     Ultra_Remember.Useed1 = Ultra_seed1;
  490.     Ultra_Remember.Useed2 = Ultra_seed2;
  491.     Ultra_Remember.Ubrw = Ultra_brw;
  492.     Ultra_Remember.Ubyt = Ultra_byt;
  493.     Ultra_Remember.Ubit = Ultra_bit;
  494.     Ultra_Remember.Uptr = Ultra_ptr;
  495.     for (i=0;i<37;i++) {
  496.         Ultra_Remember.UFSR[i] = Ultra_FSR[i];
  497.         Ultra_Remember.USWB[i] = Ultra_SWB[i];
  498.     }
  499. }
  500.  
  501. void Ultra_RecallStart()
  502. {
  503.     register    i;
  504.  
  505.     Ultra_gauss = Ultra_Remember.Ugauss;
  506.     Ultra_bits = Ultra_Remember.Ubits;
  507.     Ultra_seed1 = Ultra_Remember.Useed1;
  508.     Ultra_seed2 = Ultra_Remember.Useed2;
  509.     Ultra_brw = Ultra_Remember.Ubrw;
  510.     Ultra_byt = Ultra_Remember.Ubyt;
  511.     Ultra_bit = Ultra_Remember.Ubit;
  512.     Ultra_ptr = Ultra_Remember.Uptr;
  513.     for (i=0;i<37;i++) {
  514.         Ultra_FSR[i] = Ultra_Remember.UFSR[i];
  515.         Ultra_SWB[i] = Ultra_Remember.USWB[i];
  516.     }
  517. }
  518.  
  519. /*    Ultra_Init() computes a few global constants, initializes others and fills
  520.     in the initial Ultra_SWB array using the user-supplied seeds which must be
  521.     non-zero unsigned longs.  It terminates by calling Ultra_SaveStart() so that
  522.     you may recover the whole sequence of random numbers by calling Ultra_Init(),
  523.     with the same seeds, then calling Ultra_RecallStart(). */
  524.  
  525. Boolean Ultra_Init()
  526. {
  527.     unsigned long tempbits,ul;
  528.     short           i,j;
  529.     
  530. /* Two ulong seeds > 0 MUST be supplied by user.  There is, INTENTIONALLY, no default. */
  531.     if ((Ultra_seed1 == 0) || (Ultra_seed2 == 0)) return(false);
  532.     
  533.      for (i=0;i<37;i++) {
  534.         tempbits = 0;
  535.         for (j=32;j>0;j--) {
  536.             Ultra_seed1 *= 69069;
  537.             Ultra_seed2 ^= (Ultra_seed2 >> 15);
  538.             Ultra_seed2 ^= (Ultra_seed2 << 17);
  539.             ul = Ultra_seed1 ^ Ultra_seed2;
  540.             tempbits = (tempbits >> 1) | (0x80000000 & ul);
  541.         }
  542.         Ultra_SWB[i] = tempbits;
  543.     }
  544.     Ultra_2n31 = ((2.0/65536)/65536);
  545.     Ultra_2n63 = 0.5*Ultra_2n31*Ultra_2n31;
  546.     Ultra_2n7 = 1.0/128;
  547.     Ultra_gauss = 0.0;
  548.     Ultra_byt = Ultra_bit = 0;
  549.     Ultra_brw = ULTRABRW;            /* no borrow yet */
  550.     Ultra_SaveStart();
  551.     return(true);
  552. }
  553.